library(lattice)

ipw <- read.csv("../results/IPW.csv")
aipw <- read.csv("../results/AIPW.csv")
reg <- read.csv("../results/reg.csv")
mat <- read.csv("../results/Matching.csv")

ipw <- ipw[,2:7]
aipw <- aipw[,2:7]
reg <- reg[,2:7]
mat <- mat[,c(2,5:9)]


ipw$estim <- "IPW"
aipw$estim <- "AIPW"
reg$estim <- "Regression"
mat$estim <- "Matching"





out <- rbind(ipw, aipw, reg, mat)
out$bias <- out$ATE.hat - 5.00
out$MSE <- (out$ATE.hat - 5.00)^2

agg <- aggregate(out, by=list(nobs=out$n, alpha=out$alpha.ind,
                        beta=out$beta.ind, specif=out$spec,
                        estimat=out$estim), mean)

agg2 <- aggregate(out, by=list(nobs=out$n, alpha=out$alpha.ind,
                        beta=out$beta.ind, specif=out$spec,
                        estimat=out$estim), sd, na.rm=TRUE)

agg$SampSD <- agg2$ATE.hat


agg$RMSE <- sqrt(agg$MSE)
agg$RMSErelat <- NA
agg$RMSErelat[1:126] <- agg$RMSE[1:126] / agg$RMSE[1:126]
agg$RMSErelat[127:252] <- agg$RMSE[127:252] / agg$RMSE[1:126]
agg$RMSErelat[253:378] <- agg$RMSE[253:378] / agg$RMSE[1:126]
agg$RMSErelat[379:504] <- agg$RMSE[379:504] / agg$RMSE[1:126]



agg$specif <- as.character(agg$specif)
agg$specif[agg$specif=="1"] <- "A"
agg$specif[agg$specif=="2"] <- "B"
agg$specif[agg$specif=="3"] <- "E"
agg$specif[agg$specif=="4"] <- "F"
agg$specif[agg$specif=="5"] <- "C"
agg$specif[agg$specif=="6"] <- "G"
agg$specif[agg$specif=="7"] <- "D"

agg$alpha <- as.character(agg$alpha)
agg$alpha[agg$alpha=="1"] <- "Confounding = Low"
agg$alpha[agg$alpha=="2"] <- "Confounding = Moderate"
agg$alpha[agg$alpha=="3"] <- "Confounding = Severe"


agg$beta <- as.character(agg$beta)
agg$beta[agg$beta=="1"] <- "Outcome Mean = Linear"
agg$beta[agg$beta=="2"] <- "Outcome Mean = Nonlinear"

agg$nobs <- factor(as.character(agg$nobs), levels=c("250", "500", "1000"),
                   labels=c("n = 250", "n = 500", "n = 1000"))


mispec <- agg$specif == "E" | agg$specif == "F" | agg$specif == "G" 

trellis.device("pdf", col=FALSE)

pdf("RelativeRMSECorrect.pdf", height=11*1.5, width=8.5*1.5, bg="white")
print(dotplot(RMSErelat~specif|alpha*nobs*beta, groups=factor(estimat),
              subset=!mispec, layout=c(3,6),
              data=agg, cex=1.25,
              ylab="RMSE of Estimator for ATE Divided by RMSE of AIPW Estimator for ATE",
              xlab="Model Specification",
              col=c("white", "black", "black", "black"),
              pch=c(2,1,3,5),
              key=list(text=list(c("", "IPW", "Matching", "Regression")),
                points=list(pch=c(2,1,3,5), cex=1.25,
                  col=c("white", "black", "black", "black"))),
              panel=function(x,y,subscripts,groups,col,pch,...){
                panel.abline(h=1, col.line="gray")
                panel.dotplot(x,y,col=col[groups[subscripts]],
                              pch=pch[groups[subscripts]], ...)
              })
      )
dev.off()



mispec <- agg$specif == "E" | agg$specif == "F" | agg$specif == "G" 
pdf("BiasCorrect.pdf", height=11*1.5, width=8.5*1.5, bg="white")
print(dotplot(bias~specif|alpha*nobs*beta, groups=factor(estimat),
              subset=!mispec, layout=c(3,6),
              data=agg, cex=1.25,
              ylab="Bias of Estimator for ATE",
              xlab="Model Specification",
              col=c("black", "black", "black", "black"),
              pch=c(20,1,3,5),
              key=list(text=list(c("AIPW", "IPW", "Matching", "Regression")),
                points=list(pch=c(20,1,3,5), cex=1.25,
                  col=c("black", "black", "black", "black"))),
              panel=function(x,y,subscripts,groups,col,pch,...){
                panel.abline(h=0, col.line="gray")
                panel.dotplot(x,y,col=col[groups[subscripts]],
                              pch=pch[groups[subscripts]], ...)
              })
      )
dev.off()



mispec <- agg$specif == "E" | agg$specif == "F" | agg$specif == "G" 
pdf("RMSECorrect.pdf", height=11*1.5, width=8.5*1.5, bg="white")
print(dotplot(RMSE~specif|alpha*nobs*beta, groups=factor(estimat),
              subset=!mispec, layout=c(3,6),
              data=agg, cex=1.25,
              ylab="RMSE of Estimator for ATE",
              xlab="Model Specification",
              col=c("black", "black", "black", "black"),
              pch=c(20,1,3,5),
              key=list(text=list(c("AIPW", "IPW", "Matching", "Regression")),
                points=list(pch=c(20,1,3,5), cex=1.25,
                  col=c("black", "black", "black", "black"))),
              panel=function(x,y,subscripts,groups,col,pch,...){
                panel.abline(h=0, col.line="gray")
                panel.dotplot(x,y,col=col[groups[subscripts]],
                              pch=pch[groups[subscripts]], ...)
              })
      )
dev.off()







mispec <- agg$specif == "E" | agg$specif == "F" 
pdf("RelativeRMSEMisspec.pdf", height=11*1.5, width=8.5*1.5, bg="white")

print(dotplot(RMSErelat~specif|alpha*nobs*beta, groups=factor(estimat),
              subset=mispec, layout=c(3,6),
              data=agg, cex=1.25,
              ylab="RMSE of Estimator for ATE Divided by RMSE of AIPW Estimator for ATE",
              xlab="Model Specification",
              col=c("white", "black", "black", "black"),
              pch=c(2,1,3,5),
              key=list(text=list(c("", "IPW", "Matching", "Regression")),
                points=list(pch=c(2,1,3,5), cex=1.25,
                  col=c("white", "black", "black", "black"))),
              panel=function(x,y,subscripts,groups,col,pch,...){
                panel.abline(h=1, col.line="gray")
                panel.dotplot(x,y,col=col[groups[subscripts]],
                              pch=pch[groups[subscripts]], ...)
              })
      )
dev.off()

mispec <- agg$specif == "E" | agg$specif == "F" 

pdf("BiasMisspec.pdf", height=11*1.5, width=8.5*1.5, bg="white")
print(dotplot(bias~specif|alpha*nobs*beta, groups=factor(estimat),
              subset=mispec, layout=c(3,6),
              data=agg, cex=1.25,
              ylab="Bias of Estimator for ATE",
              xlab="Model Specification",
              col=c("black", "black", "black", "black"),
              pch=c(20,1,3,5),
              key=list(text=list(c("AIPW", "IPW", "Matching", "Regression")),
                points=list(pch=c(20,1,3,5), cex=1.25,
                  col=c("black", "black", "black", "black"))),
              panel=function(x,y,subscripts,groups,col,pch,...){
                panel.abline(h=0, col.line="gray")
                panel.dotplot(x,y,col=col[groups[subscripts]],
                              pch=pch[groups[subscripts]], ...)
              })
      )
dev.off()


pdf("RMSEMisspec.pdf", height=11*1.5, width=8.5*1.5, bg="white")
print(dotplot(RMSE~specif|alpha*nobs*beta, groups=factor(estimat),
              subset=mispec, layout=c(3,6),
              data=agg, cex=1.25,
              ylab="RMSE of Estimator for ATE",
              xlab="Model Specification",
              col=c("black", "black", "black", "black"),
              pch=c(20,1,3,5),
              key=list(text=list(c("AIPW", "IPW", "Matching", "Regression")),
                points=list(pch=c(20,1,3,5), cex=1.25,
                  col=c("black", "black", "black", "black"))),
              panel=function(x,y,subscripts,groups,col,pch,...){
                panel.abline(h=0, col.line="gray")
                panel.dotplot(x,y,col=col[groups[subscripts]],
                              pch=pch[groups[subscripts]], ...)
              })
      )
dev.off()











